/*
 * Voronoi.h
 *
 *  Created on: Jun 20, 2011
 *      Author: marchi
 */

#ifndef VORONOI_H_
#define VORONOI_H_
#include <string>
#include <vector>
#include <map>
#include <list>
#include "FAtoms.h"
#include "ResidueCM.h"
#include "ResidueTypes.h"
#include <cstring>
#include "Array.h"
#include <fstream>
#include <iomanip>
#include "VoronoiPrint.h"
#include "Topol.h"
#include "isTypes.h"

using namespace Array;

using namespace std;
#include "voro++.hh"
using namespace voro;

using namespace std;
const int NNN=26;
const char h[]="Hxxxxx";
const double SMALL=1.0e-4;


// Golden ratio constants
const double Phi=0.5*(1+sqrt(5.0));
const double phi=0.5*(1-sqrt(5.0));

// Set up the number of blocks that the container is divided
// into.
const int n_x=5,n_y=5,n_z=5;


class wall_initial_shape : public wall {
	void putR0(double rsp){
		// Initialize the Voronoi cell to be a cube of side length 2, centered
	     // on the origin

		voronoicell_neighbor v;
		double a=2.0*rsp;
		double myrsq=a*a;
		v.init(-2*a,2*a,-2*a,2*a,-2*a,2*a);

		// Cut the cell by 250 random planes which are all a distance a away
		// from the origin, to make an approximation to a sphere
		for(int i=0;i<250;i++) {
			double x=2*double(rand())/RAND_MAX-1;
			double y=2*double(rand())/RAND_MAX-1;
			double z=2*double(rand())/RAND_MAX-1;
			double rsq=x*x+y*y+z*z;
			double r=1/sqrt(rsq);x*=r;y*=r;z*=r;
			x*=a;y*=a;z*=a;
			v.nplane(x,y,z,myrsq,-1);
		}
		vor[rsp]=new voronoicell_neighbor;
		*vor[rsp]=v;
		cout << " Wall No. " << count << " defined of radius " << rsp  <<endl;
		count++;
	}
public:
	wall_initial_shape(): R0(0.0){};
	wall_initial_shape(const double rsp){
		R0=rsp;
	}
	void operator()(double rsp){
		putR0(rsp);
	};
	void setR0(double R){
		if(vor.count(R) > 0) {
			R0=R;
		} else{
			putR0(R);
			R0=R;
		}
	};
	bool point_inside(double x,double y,double z) {return true;}
	bool cut_cell(voronoicell &c,double x,double y,double z) {

			// Set the cell to be equal to the dodecahedron
		return true;
	}
	bool cut_cell(voronoicell_neighbor &c,double x,double y,double z) {

			// Set the cell to be equal to the dodecahedron
		c=*vor[R0];
		return true;
	}
	virtual ~wall_initial_shape(){
		for(map<double,voronoicell_neighbor *>::iterator it=vor.begin();it != vor.end();it++)
			delete it->second;
		vor.clear();
		count=0;
	}

private:
	map<double,voronoicell_neighbor *> vor;
		double R0;
		static int count;
};
class Voronoi {
protected:
	matrix co;
	matrix oc;
	static bool Spheres;
	static int nresid,nr,nc;
	static float time;
	static int * types;
	string * TypesName;
	double * Vol;
	array1<double> Vols;
	vector<Dvect> Centroids;
	vector<Dvect> Cxcm;
	array2<double> area;
	static container_periodic_poly * MyconP;
	static container_poly * MyconNP;
	static wall_initial_shape * wis0;

	particle_order * porder;
	static string * label;
	vector<int> cindex;
	vector<int> * Neighs;
	vector<double> * Surface;
	vector<int> RealResidue;
	vector<vector<double> > Phobic,Philic,Mixt;
	vector<string> unique_domains;
	void gather(vector<int> & y);
	virtual void WriteIt(std::ofstream &);
	vector<int> ResType;
	vector<string> restype_s;
	vector<int> restype_i;
	Voronoi();
public:
	Voronoi(int, char ***,bool);
	Voronoi(int, const vector<string> &,bool);
	void StartP(float, Atoms &);
	void StartNP(float, Atoms &);
	void StartNP(Atoms &);
	void getData(vector<Topol_NS::MySubunit> &);
	virtual void getDelaunay(Cavities::isTypes &);
	virtual void sTname(Topol_NS::Topol & y ){};
	virtual vector<double> & gTotHG(){static vector<double> t;return t;};
	virtual void sTotHG(vector<double> & y){}
	virtual void copyhgtotot(){};
	virtual void copytottohg(){};
	virtual void Channels(){};

	const matrix & getCO(){return co;}
	const matrix & getOC(){return oc;}
	const array1<double> & getVols(){return Vols;}
	const vector<Dvect> & getCentroids(){return Centroids;}
	const vector<int> & getCindex(){return cindex;}
	void testVol();
	string & name(int i) const {return label[i];}
	void setTypes(int,const AtomIndex [],const string *);
	double getVolN(int n){return Vols[n];};
	int getnresid(){return nresid;}
	int getnr(){return nr;}
	void testit();
	void setResTypes(Topol_NS::Topol &);
	void setSpheres(){Spheres=true;}
	int getTypes(int i) {return types[i];};
	int getCentroidsN(){return Centroids.size();}
	string & getTypesName(int i) const {return TypesName[types[i]];};

	int getTypesRes(int nc) {
		int First=ResidueCM::getind(nc)[0];
		return getTypes(First);
	}
	int VoronoiCellID(const double & xc,const double & yc,const double & zc){
		double rx,ry,rz;
		int iD;
		cout << " Here " <<endl;
		if(MyconP) cout << " MyconP  "<< endl;
		if(MyconNP) cout << " MyconNP  "<< endl;
		exit(1);
		try{if(!MyconP->find_voronoi_cell(xc,yc,zc,rx,ry,rz,iD))
			throw "Error while finding Voronoi cell ";}
		catch(const char * s) {std::cout << s << std::endl;exit(1);}
		return cindex[iD];
	}
	int VoronoiCellID(const double & xc,const double & yc,const double & zc,
			double & rx, double & ry, double & rz){
		int iD;
		cout << " Here 2" <<endl;
		if(MyconP) cout << " MyconP  "<< endl;
		if(MyconNP) cout << " MyconNP  "<< endl;
		exit(1);
		if(MyconP){
			try{if(!MyconP->find_voronoi_cell(xc,yc,zc,rx,ry,rz,iD))
				throw "Error while finding Voronoi cell ";}
			catch(const char * s) {std::cout << s << std::endl;exit(1);}
		} else if(MyconNP) {
			try{if(!MyconNP->find_voronoi_cell(xc,yc,zc,rx,ry,rz,iD))
				throw "Error while finding Voronoi cell ";}
			catch(const char * s) {std::cout << s << std::endl;exit(1);}
		}
		return cindex[iD];
	}
	virtual ~Voronoi(){
		cout << " delete voronoi " <<endl;
			delete [] TypesName;
			delete [] Vol;
			delete [] Neighs;
			delete [] Surface;
			delete MyconP;
			delete MyconNP;
			delete porder;
		};
	friend std::ofstream & operator<<(std::ofstream & fout, Voronoi  & y){
		y.WriteIt(fout);
		return fout;
	}
};
#endif /* VORONOI_H_ */
